Analysis of GLDS-208_Array from NASA GeneLab

This R markdown file was auto-generated by the iDEP website Using iDEP 0.91, originally by Steven

Ge SX, Son EW, Yao R: iDEP: an integrated web application for differential expression and pathway analysis of RNA-Seq data. BMC Bioinformatics 2018, 19(1):534. PMID:30567491

1. Read data

First we set up the working directory to where the files are saved.

 setwd('~/Documents/HTML_R/GLDS208_Array')

R packages and iDEP core Functions. Users can also download the iDEP_core_functions.R file. Many R packages needs to be installed first. This may take hours. Each of these packages took years to develop.So be a patient thief. Sometimes dependencies needs to be installed manually. If you are using an older version of R, and having trouble with package installation, try un-install the current version of R, delete all folders and files (C:/Program Files/R/R-3.4.3), and reinstall from scratch.

 if(file.exists('iDEP_core_functions.R'))
    source('iDEP_core_functions.R') else 
    source('https://raw.githubusercontent.com/iDEP-SDSU/idep/master/shinyapps/idep/iDEP_core_functions.R')

We are using the downloaded gene expression file where gene IDs has been converted to Ensembl gene IDs. This is because the ID conversion database is too large to download. You can use your original file if your file uses Ensembl ID, or you do not want to use the pathway files available in iDEP (or it is not available).

 inputFile <- 'GLDS208_Array_Expression.csv'
 sampleInfoFile <- 'GLDS208_Array_Sampleinfo.csv'
 gldsMetadataFile <- 'GLDS208_Array_Metadata.csv'
 geneInfoFile <- 'Arabidopsis_thaliana__athaliana_eg_gene_GeneInfo.csv' #Gene symbols, location etc. 
 geneSetFile <- 'Arabidopsis_thaliana__athaliana_eg_gene.db'  # pathway database in SQL; can be GMT format 
 STRING10_speciesFile <- 'https://raw.githubusercontent.com/iDEP-SDSU/idep/master/shinyapps/idep/STRING10_species.csv'

Parameters for reading data

 input_missingValue <- 'geneMedian' #Missing values imputation method
 input_dataFileFormat <- 1  #1- read counts, 2 FKPM/RPKM or DNA microarray
 input_minCounts <- 0.5 #Min counts
 input_NminSamples <- 1 #Minimum number of samples 
 input_countsLogStart <- 4  #Pseudo count for log CPM
 input_CountsTransform <- 1 #Methods for data transformation of counts. 1-EdgeR's logCPM 2-VST, 3-rlog 
readMetadata.out <- readMetadata(gldsMetadataFile)
library(knitr)   #  install if needed. for showing tables with kable
library(kableExtra)
kable( readMetadata.out ) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
ZoneI_Rep1 ZoneI_Rep2 ZoneI_Rep3 ZoneII_Rep1 ZoneII_Rep2 ZoneII_Rep3
Sample.LongId Atha.Col.0.Root.Zone.I.Rep1.Array Atha.Col.0.Root.Zone.I.Rep2.Array Atha.Col.0.Root.Zone.I.Rep3.Array Atha.Col.0.Root.Zone.II.Rep1.Array Atha.Col.0.Root.Zone.II.Rep2.Array Atha.Col.0.Root.Zone.II.Rep3.Array
Sample.Id Atha.Col.0.Root.Zone.I.Rep1 Atha.Col.0.Root.Zone.I.Rep2 Atha.Col.0.Root.Zone.I.Rep3 Atha.Col.0.Root.Zone.II.Rep1 Atha.Col.0.Root.Zone.II.Rep2 Atha.Col.0.Root.Zone.II.Rep3
Sample.Name Atha_Col-0_Root_Zone-I_Rep1 Atha_Col-0_Root_Zone-I_Rep2 Atha_Col-0_Root_Zone-I_Rep3 Atha_Col-0_Root_Zone-II_Rep1 Atha_Col-0_Root_Zone-II_Rep2 Atha_Col-0_Root_Zone-II_Rep3
GLDS 208 208 208 208 208 208
Accession GLDS-208 GLDS-208 GLDS-208 GLDS-208 GLDS-208 GLDS-208
Hardware Petri dish Petri dish Petri dish Petri dish Petri dish Petri dish
Tissue Meristem zone 1 Meristem zone 1 Meristem zone 1 Elongation zone 2 Elongation zone 2 Elongation zone 2
Age 8 days 8 days 8 days 8 days 8 days 8 days
Organism Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana
Ecotype Col-0 Col-0 Col-0 Col-0 Col-0 Col-0
Genotype WT WT WT WT WT WT
Variety Col-0 WT Col-0 WT Col-0 WT Col-0 WT Col-0 WT Col-0 WT
Radiation Background Earth Background Earth Background Earth Background Earth Background Earth Background Earth
Gravity Terrestrial Terrestrial Terrestrial Terrestrial Terrestrial Terrestrial
Developmental 8 day old seedlings root tip meristem and elongation zone 8 day old seedlings root tip meristem and elongation zone 8 day old seedlings root tip meristem and elongation zone 8 day old seedlings root tip meristem and elongation zone 8 day old seedlings root tip meristem and elongation zone 8 day old seedlings root tip meristem and elongation zone
Time.series.or.Concentration.gradient Single time point Single time point Single time point Single time point Single time point Single time point
Light White light White light White light White light White light White light
Assay..RNAseq. RNAseq Transcription Profiling RNAseq Transcription Profiling RNAseq Transcription Profiling RNAseq Transcription Profiling RNAseq Transcription Profiling RNAseq Transcription Profiling
Temperature 18-20 18-20 18-20 18-20 18-20 18-20
Treatment.type Comparative gene expression analysis in the Arabidopsis thaliana root apex using RNA-seq and microarray transcriptome profiles (microarray) Comparative gene expression analysis in the Arabidopsis thaliana root apex using RNA-seq and microarray transcriptome profiles (microarray) Comparative gene expression analysis in the Arabidopsis thaliana root apex using RNA-seq and microarray transcriptome profiles (microarray) Comparative gene expression analysis in the Arabidopsis thaliana root apex using RNA-seq and microarray transcriptome profiles (microarray) Comparative gene expression analysis in the Arabidopsis thaliana root apex using RNA-seq and microarray transcriptome profiles (microarray) Comparative gene expression analysis in the Arabidopsis thaliana root apex using RNA-seq and microarray transcriptome profiles (microarray)
Treatment.intensity x x x x x x
Treament.timing x x x x x x 
 readData.out <- readData(inputFile) 
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
   kable( head(readData.out$data) ) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
ZoneI_Rep1 ZoneI_Rep2 ZoneI_Rep3 ZoneII_Rep1 ZoneII_Rep2 ZoneII_Rep3
AT2G36100 2.807355 2.584963 2.807355 3.906891 3.906891 3.807355
AT5G66390 2.584963 2.584963 2.584963 3.906891 3.807355 3.807355
AT5G04960 3.169925 2.584963 2.584963 3.906891 3.906891 3.906891
AT5G67400 2.807355 2.807355 2.584963 3.906891 3.906891 3.906891
AT5G42180 2.584963 2.584963 2.584963 3.807355 3.807355 3.700440
AT5G35190 2.584963 2.584963 2.321928 3.700440 3.700440 3.700440 
 readSampleInfo.out <- readSampleInfo(sampleInfoFile) 
 kable( readSampleInfo.out ) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
Tissue
ZoneI_Rep1 Meristem_zone1
ZoneI_Rep2 Meristem_zone1
ZoneI_Rep3 Meristem_zone1
ZoneII_Rep1 Elongation_zone2
ZoneII_Rep2 Elongation_zone2
ZoneII_Rep3 Elongation_zone2 
 input_selectOrg ="NEW" 
 input_selectGO <- 'GOBP'   #Gene set category 
 input_noIDConversion = TRUE  
 allGeneInfo.out <- geneInfo(geneInfoFile) 
 converted.out = NULL 
 convertedData.out <- convertedData()    
 nGenesFilter()  
## [1] "16156 genes in 6 samples. 16156  genes passed filter.\n Original gene IDs used."
 convertedCounts.out <- convertedCounts()  # converted counts, just for compatibility

2. Pre-process

# Read counts per library 
 parDefault = par() 
 par(mar=c(12,4,2,2)) 
 # barplot of total read counts
 x <- readData.out$rawCounts
 groups = as.factor( detectGroups(colnames(x ) ) )
 if(nlevels(groups)<=1 | nlevels(groups) >20 )  
  col1 = 'green'  else
  col1 = rainbow(nlevels(groups))[ groups ]             
         
 barplot( colSums(x)/1e6, 
        col=col1,las=3, main="Total read counts (millions)")  

 readCountsBias()  # detecting bias in sequencing depth 
## [1] 0.5852447
## [1] 0.5852447
## [1] "No bias detected"
 # Box plot 
 x = readData.out$data 
 boxplot(x, las = 2, col=col1,
    ylab='Transformed expression levels',
    main='Distribution of transformed data') 

 #Density plot 
 par(parDefault) 
## Warning in par(parDefault): graphical parameter "cin" cannot be set
## Warning in par(parDefault): graphical parameter "cra" cannot be set
## Warning in par(parDefault): graphical parameter "csi" cannot be set
## Warning in par(parDefault): graphical parameter "cxy" cannot be set
## Warning in par(parDefault): graphical parameter "din" cannot be set
## Warning in par(parDefault): graphical parameter "page" cannot be set
 densityPlot()       

 # Scatter plot of the first two samples 
 plot(x[,1:2],xlab=colnames(x)[1],ylab=colnames(x)[2], 
    main='Scatter plot of first two samples') 

 ####plot gene or gene family
 input_selectOrg ="BestMatch" 
 input_geneSearch <- 'HOXA' #Gene ID for searching 
 genePlot()  
## NULL
 input_useSD <- 'FALSE' #Use standard deviation instead of standard error in error bar? 
 geneBarPlotError()       
## NULL

3. Heatmap

 # hierarchical clustering tree
 x <- readData.out$data
 maxGene <- apply(x,1,max)
 # remove bottom 25% lowly expressed genes, which inflate the PPC
 x <- x[which(maxGene > quantile(maxGene)[1] ) ,] 
 plot(as.dendrogram(hclust2( dist2(t(x)))), ylab="1 - Pearson C.C.", type = "rectangle") 

 #Correlation matrix
 input_labelPCC <- TRUE #Show correlation coefficient? 
 correlationMatrix() 

 # Parameters for heatmap
 input_nGenes <- 1000   #Top genes for heatmap
 input_geneCentering <- TRUE    #centering genes ?
 input_sampleCentering <- FALSE #Center by sample?
 input_geneNormalize <- FALSE   #Normalize by gene?
 input_sampleNormalize <- FALSE #Normalize by sample?
 input_noSampleClustering <- FALSE  #Use original sample order
 input_heatmapCutoff <- 4   #Remove outliers beyond number of SDs 
 input_distFunctions <- 1   #which distant funciton to use
 input_hclustFunctions <- 1 #Linkage type
 input_heatColors1 <- 1 #Colors
 input_selectFactorsHeatmap <- 'Tissue' #Sample coloring factors 
 png('heatmap.png', width = 10, height = 15, units = 'in', res = 300) 
 staticHeatmap() 
 dev.off()  
## png 
##   2

[heatmap] (heatmap.png)

 heatmapPlotly() # interactive heatmap using Plotly

4. K-means clustering

 input_nGenesKNN <- 2000    #Number of genes fro k-Means
 input_nClusters <- 4   #Number of clusters 
 maxGeneClustering = 12000
 input_kmeansNormalization <- 'geneMean'    #Normalization
 input_KmeansReRun <- 0 #Random seed 

 distributionSD()  #Distribution of standard deviations 

 KmeansNclusters()  #Number of clusters 

 Kmeans.out = Kmeans()   #Running K-means 
 KmeansHeatmap()   #Heatmap for k-Means 

 #Read gene sets for enrichment analysis 
 sqlite  <- dbDriver('SQLite')
 input_selectGO3 <- 'GOBP'  #Gene set category
 input_minSetSize <- 15 #Min gene set size
 input_maxSetSize <- 2000   #Max gene set size 
 GeneSets.out <-readGeneSets( geneSetFile,
    convertedData.out, input_selectGO3,input_selectOrg,
    c(input_minSetSize, input_maxSetSize)  )  
 # Alternatively, users can use their own GMT files by
 #GeneSets.out <- readGMTRobust('somefile.GMT')  
 results <- KmeansGO()  #Enrichment analysis for k-Means clusters   
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
Cluster adj.Pval Genes Pathways
A 8.85e-20 125 Response to organic substance
1.86e-17 106 Response to hormone
1.86e-17 99 Oxidation-reduction process
4.90e-17 106 Response to endogenous stimulus
5.84e-17 123 Response to abiotic stimulus
2.24e-16 101 Response to oxygen-containing compound
1.24e-15 96 Cellular response to chemical stimulus
2.66e-15 82 Response to acid chemical
1.40e-14 87 Response to external stimulus
7.95e-14 55 Cell wall organization or biogenesis
B 5.85e-25 48 Cell wall organization or biogenesis
3.43e-22 41 External encapsulating structure organization
3.43e-22 40 Cell wall organization
3.07e-18 25 Drug catabolic process
1.88e-17 20 Phenylpropanoid metabolic process
9.55e-16 22 Cellular response to toxic substance
7.00e-15 25 Plant-type cell wall organization or biogenesis
5.01e-14 15 Phenylpropanoid biosynthetic process
1.36e-13 19 Cellular oxidant detoxification
1.48e-13 15 Hydrogen peroxide catabolic process
C 6.78e-20 38 Cell cycle
6.78e-20 31 Cell division
4.74e-19 26 Mitotic cell cycle process
1.04e-18 28 Mitotic cell cycle
1.05e-17 30 Cell cycle process
3.20e-15 19 Regulation of mitotic cell cycle
2.08e-13 22 Regulation of cell cycle
3.33e-13 17 Regulation of cell cycle process
3.36e-12 15 Cell proliferation
4.83e-12 18 Nuclear division
D 3.70e-82 96 Ribosome biogenesis
1.86e-76 100 Ribonucleoprotein complex biogenesis
1.96e-74 80 RRNA metabolic process
3.80e-74 91 NcRNA processing
5.27e-74 79 RRNA processing
1.11e-70 94 NcRNA metabolic process
1.03e-55 104 RNA processing
2.38e-55 125 Cellular component biogenesis
1.24e-31 34 Ribosomal large subunit biogenesis
1.11e-21 20 Maturation of LSU-rRNA 
 input_seedTSNE <- 0    #Random seed for t-SNE
 input_colorGenes <- TRUE   #Color genes in t-SNE plot? 
 tSNEgenePlot()  #Plot genes using t-SNE

5. PCA and beyond

 input_selectFactors <- 'Sample_Name'   #Factor coded by color
 input_selectFactors2 <- 'Sample_Name'  #Factor coded by shape
 input_tsneSeed2 <- 0   #Random seed for t-SNE 
 #PCA, MDS and t-SNE plots
 PCAplot()  

 MDSplot() 

 tSNEplot()  

 #Read gene sets for pathway analysis using PGSEA on principal components 
 input_selectGO6 <- 'GOBP' 
 GeneSets.out <-readGeneSets( geneSetFile,
    convertedData.out, input_selectGO6,input_selectOrg,
    c(input_minSetSize, input_maxSetSize)  )  
 PCApathway() # Run PGSEA analysis 
## Warning: Package 'KEGG.db' is deprecated and will be removed from Bioconductor
##   version 3.12

 cat( PCA2factor() )   #The correlation between PCs with factors 
## 
##  Correlation between Principal Components (PCs) with factors
## PC1 is correlated with Tissue (p=7.75e-05).

6. DEG1

 input_CountsDEGMethod <- 2 #DESeq2= 3,limma-voom=2,limma-trend=1 
 input_limmaPval <- 0.1 #FDR cutoff
 input_limmaFC <- 2 #Fold-change cutoff
 input_selectModelComprions <- 'Tissue: Elongation_zone2 vs. Meristem_zone1'    #Selected comparisons
 input_selectFactorsModel <- 'Tissue'   #Selected comparisons
 input_selectInteractions <- NULL   #Selected comparisons
 input_selectBlockFactorsModel <- NULL  #Selected comparisons
 factorReferenceLevels.out <- c('Tissue:Meristem_zone1') 

 limma.out <- limma()
 DEG.data.out <- DEG.data()
 limma.out$comparisons 
## [1] "ZoneII-ZoneI"
 input_selectComparisonsVenn = limma.out$comparisons[1:3] # use first three comparisons
 input_UpDownRegulated <- FALSE #Split up and down regulated genes 
 vennPlot() # Venn diagram 

  sigGeneStats() # number of DEGs as figure 

  sigGeneStatsTable() # number of DEGs as table 
##               Comparisons  Up Down
## ZoneII-ZoneI ZoneII-ZoneI 538  126

7. DEG2

 input_selectContrast = limma.out$comparisons[1] # use first  comparisons 
 selectedHeatmap.data.out <- selectedHeatmap.data()
## Error in findContrastSamples(input_selectContrast, colnames(convertedData.out), : object 'c.out' not found
 selectedHeatmap()   # heatmap for DEGs in selected comparison
## Error in selectedHeatmap(): object 'selectedHeatmap.data.out' not found
 # Save gene lists and data into files
 write.csv( selectedHeatmap.data()$genes, 'heatmap.data.csv') 
## Error in findContrastSamples(input_selectContrast, colnames(convertedData.out), : object 'c.out' not found
 write.csv(DEG.data(),'DEG.data.csv' )
 write(AllGeneListsGMT() ,'AllGeneListsGMT.gmt')
 input_selectGO2 <- 'GOBP'  #Gene set category 
 geneListData.out <- geneListData()  
 volcanoPlot()  

  scatterPlot()  
## Error in findContrastSamples(input_selectContrast, colnames(convertedData.out), : object 'c.out' not found
  MAplot()  
## Error in findContrastSamples(input_selectContrast, colnames(convertedData.out), : object 'c.out' not found
  geneListGOTable.out <- geneListGOTable()  
## Error in geneListGOTable(): object 'selectedHeatmap.data.out' not found
 # Read pathway data again 
 GeneSets.out <-readGeneSets( geneSetFile,
    convertedData.out, input_selectGO2,input_selectOrg,
    c(input_minSetSize, input_maxSetSize)  ) 
 input_removeRedudantSets <- TRUE   #Remove highly redundant gene sets? 
 results <- geneListGO()  #Enrichment analysis
## Error in geneListGO(): object 'geneListGOTable.out' not found
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
Cluster adj.Pval Genes Pathways
A 8.85e-20 125 Response to organic substance
1.86e-17 106 Response to hormone
1.86e-17 99 Oxidation-reduction process
4.90e-17 106 Response to endogenous stimulus
5.84e-17 123 Response to abiotic stimulus
2.24e-16 101 Response to oxygen-containing compound
1.24e-15 96 Cellular response to chemical stimulus
2.66e-15 82 Response to acid chemical
1.40e-14 87 Response to external stimulus
7.95e-14 55 Cell wall organization or biogenesis
B 5.85e-25 48 Cell wall organization or biogenesis
3.43e-22 41 External encapsulating structure organization
3.43e-22 40 Cell wall organization
3.07e-18 25 Drug catabolic process
1.88e-17 20 Phenylpropanoid metabolic process
9.55e-16 22 Cellular response to toxic substance
7.00e-15 25 Plant-type cell wall organization or biogenesis
5.01e-14 15 Phenylpropanoid biosynthetic process
1.36e-13 19 Cellular oxidant detoxification
1.48e-13 15 Hydrogen peroxide catabolic process
C 6.78e-20 38 Cell cycle
6.78e-20 31 Cell division
4.74e-19 26 Mitotic cell cycle process
1.04e-18 28 Mitotic cell cycle
1.05e-17 30 Cell cycle process
3.20e-15 19 Regulation of mitotic cell cycle
2.08e-13 22 Regulation of cell cycle
3.33e-13 17 Regulation of cell cycle process
3.36e-12 15 Cell proliferation
4.83e-12 18 Nuclear division
D 3.70e-82 96 Ribosome biogenesis
1.86e-76 100 Ribonucleoprotein complex biogenesis
1.96e-74 80 RRNA metabolic process
3.80e-74 91 NcRNA processing
5.27e-74 79 RRNA processing
1.11e-70 94 NcRNA metabolic process
1.03e-55 104 RNA processing
2.38e-55 125 Cellular component biogenesis
1.24e-31 34 Ribosomal large subunit biogenesis
1.11e-21 20 Maturation of LSU-rRNA

STRING-db API access. We need to find the taxonomy id of your species, this used by STRING. First we try to guess the ID based on iDEP’s database. Users can also skip this step and assign NCBI taxonomy id directly by findTaxonomyID.out = 10090 # mouse 10090, human 9606 etc.

 STRING10_species = read.csv(STRING10_speciesFile)  
 ix = grep('Arabidopsis thaliana', STRING10_species$official_name ) 
 findTaxonomyID.out <- STRING10_species[ix,1] # find taxonomyID
 findTaxonomyID.out  
## [1] 3702

Enrichment analysis using STRING

  STRINGdb_geneList.out <- STRINGdb_geneList() #convert gene lists
 input_STRINGdbGO <- 'Process'  #'Process', 'Component', 'Function', 'KEGG', 'Pfam', 'InterPro' 
 results <- stringDB_GO_enrichmentData()  # enrichment using STRING 
## Error in stringDB_GO_enrichmentData(): object 'selectedHeatmap.data.out' not found
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
Cluster adj.Pval Genes Pathways
A 8.85e-20 125 Response to organic substance
1.86e-17 106 Response to hormone
1.86e-17 99 Oxidation-reduction process
4.90e-17 106 Response to endogenous stimulus
5.84e-17 123 Response to abiotic stimulus
2.24e-16 101 Response to oxygen-containing compound
1.24e-15 96 Cellular response to chemical stimulus
2.66e-15 82 Response to acid chemical
1.40e-14 87 Response to external stimulus
7.95e-14 55 Cell wall organization or biogenesis
B 5.85e-25 48 Cell wall organization or biogenesis
3.43e-22 41 External encapsulating structure organization
3.43e-22 40 Cell wall organization
3.07e-18 25 Drug catabolic process
1.88e-17 20 Phenylpropanoid metabolic process
9.55e-16 22 Cellular response to toxic substance
7.00e-15 25 Plant-type cell wall organization or biogenesis
5.01e-14 15 Phenylpropanoid biosynthetic process
1.36e-13 19 Cellular oxidant detoxification
1.48e-13 15 Hydrogen peroxide catabolic process
C 6.78e-20 38 Cell cycle
6.78e-20 31 Cell division
4.74e-19 26 Mitotic cell cycle process
1.04e-18 28 Mitotic cell cycle
1.05e-17 30 Cell cycle process
3.20e-15 19 Regulation of mitotic cell cycle
2.08e-13 22 Regulation of cell cycle
3.33e-13 17 Regulation of cell cycle process
3.36e-12 15 Cell proliferation
4.83e-12 18 Nuclear division
D 3.70e-82 96 Ribosome biogenesis
1.86e-76 100 Ribonucleoprotein complex biogenesis
1.96e-74 80 RRNA metabolic process
3.80e-74 91 NcRNA processing
5.27e-74 79 RRNA processing
1.11e-70 94 NcRNA metabolic process
1.03e-55 104 RNA processing
2.38e-55 125 Cellular component biogenesis
1.24e-31 34 Ribosomal large subunit biogenesis
1.11e-21 20 Maturation of LSU-rRNA

PPI network retrieval and analysis

 input_nGenesPPI <- 100 #Number of top genes for PPI retrieval and analysis 
 stringDB_network1(1) #Show PPI network

Generating interactive PPI

 write(stringDB_network_link(), 'PPI_results.html') # write results to html file 
## Warning: 'string_db$get_link' is deprecated.
## Use 'Contact developers to request functionality' instead.
## See help("Deprecated")

## Warning: 'string_db$get_link' is deprecated.
## Use 'Contact developers to request functionality' instead.
## See help("Deprecated")

## Warning: 'string_db$get_link' is deprecated.
## Use 'Contact developers to request functionality' instead.
## See help("Deprecated")
 browseURL('PPI_results.html') # open in browser

8. Pathway analysis

 input_selectContrast1 = limma.out$comparisons[1] 
 #input_selectContrast1 = limma.out$comparisons[3] # manually set
 input_selectGO <- 'GOBP'   #Gene set category 
 #input_selectGO='custom' # if custom gmt file
 input_minSetSize <- 15 #Min size for gene set
 input_maxSetSize <- 2000   #Max size for gene set 
 # Read pathway data again 
 GeneSets.out <-readGeneSets( geneSetFile,
    convertedData.out, input_selectGO,input_selectOrg,
    c(input_minSetSize, input_maxSetSize)  ) 
 input_pathwayPvalCutoff <- 0.2 #FDR cutoff
 input_nPathwayShow <- 30   #Top pathways to show
 input_absoluteFold <- FALSE    #Use absolute values of fold-change?
 input_GenePvalCutoff <- 1  #FDR to remove genes 

 input_pathwayMethod = 1  # 1  GAGE
 gagePathwayData.out <- gagePathwayData()  # pathway analysis using GAGE  
   
 results <- gagePathwayData.out  #Enrichment analysis for k-Means clusters  
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
Direction GAGE analysis: ZoneII vs ZoneI statistic Genes adj.Pval
Down Ribonucleoprotein complex biogenesis -15.5111 438 1.7e-43
NcRNA metabolic process -15.0653 425 7.0e-42
Ribosome biogenesis -14.9711 343 6.9e-40
NcRNA processing -14.6745 357 2.7e-39
RRNA processing -12.9994 239 1.0e-29
RRNA metabolic process -12.9212 244 1.1e-29
RNA modification -11.0436 321 2.0e-23
Macromolecule methylation -9.9078 211 1.5e-18
Chromatin organization -9.1348 384 8.4e-17
MRNA metabolic process -9.0162 481 1.4e-16
MRNA processing -8.9205 378 5.8e-16
TRNA metabolic process -8.8598 170 6.3e-15
Ribosomal large subunit biogenesis -8.4945 99 1.5e-12
Cell cycle process -8.4921 396 7.3e-15
RNA splicing -8.2441 265 1.7e-13
Methylation -8.1162 364 1.4e-13
Ribonucleoprotein complex assembly -8.0979 174 1.5e-12
Ribonucleoprotein complex subunit organization -8.0957 182 1.5e-12
Nuclear division -7.8255 189 2.5e-12
Cellular response to DNA damage stimulus -7.8025 337 1.5e-12
Organelle fission -7.7706 226 2.5e-12
DNA repair -7.7373 314 2.3e-12
TRNA processing -7.5733 108 6.2e-11
Ribosomal small subunit biogenesis -7.2565 80 1.2e-09
RNA splicing, via transesterification reactions -7.1742 194 2.5e-10
RNA splicing, via transesterification reactions with bulged adenosine as nucleophile -7.1742 194 2.5e-10
DNA replication -7.1335 162 2.5e-10
Mitotic cell cycle process -7.0912 235 1.9e-10
Mitotic cell cycle phase transition -6.9779 76 3.0e-09
DNA conformation change -6.8825 142 1.5e-09 
 pathwayListData.out = pathwayListData() 
 enrichmentPlot(pathwayListData.out, 25  ) 

  enrichmentNetwork(pathwayListData.out )  

  enrichmentNetworkPlotly(pathwayListData.out) 

## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
 input_pathwayMethod = 3  # 1  fgsea 
 fgseaPathwayData.out <- fgseaPathwayData() #Pathway analysis using fgsea 
## Warning in fgsea(pathways = gmt, stats = fold, minSize = input_minSetSize, :
## You are trying to run fgseaSimple. It is recommended to use fgseaMultilevel. To
## run fgseaMultilevel, you need to remove the nperm argument in the fgsea function
## call.
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (84.71% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaSimple(...): There were 1 pathways for which P-values were not
## calculated properly due to unbalanced gene-level statistic values
 results <- fgseaPathwayData.out  #Enrichment analysis for k-Means clusters 
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
Direction GSEA analysis: ZoneII vs ZoneI NES Genes adj.Pval
Down Ribosome biogenesis -3.4773 343 1.2e-02
NcRNA processing -3.4292 357 1.4e-02
Ribonucleoprotein complex biogenesis -3.4005 438 2.1e-02
RRNA metabolic process -3.3857 244 7.1e-03
RRNA processing -3.3849 239 7.1e-03
NcRNA metabolic process -3.3779 425 2.0e-02
RNA processing -3.1844 780 1.4e-01
RNA modification -3.1735 321 1.1e-02
Ribosomal large subunit biogenesis -3.0694 99 3.6e-03
Macromolecule methylation -3.0673 211 6.0e-03
Mitotic cell cycle phase transition -2.9952 76 3.4e-03
Ribonucleoprotein complex assembly -2.9439 174 5.0e-03
RNA methylation -2.9384 64 3.4e-03
Ribonucleoprotein complex subunit organization -2.9247 182 5.3e-03
Cell cycle phase transition -2.9146 81 3.4e-03
TRNA metabolic process -2.9091 170 4.8e-03
Chromosome organization -2.881 569 4.3e-02
TRNA processing -2.8631 108 3.8e-03
Ribosomal small subunit biogenesis -2.8217 80 3.4e-03
TRNA modification -2.8189 78 3.4e-03
Regulation of mitotic nuclear division -2.8038 63 3.4e-03
Nuclear division -2.7958 189 5.4e-03
Regulation of cell cycle process -2.7932 136 4.4e-03
Mitochondrial RNA metabolic process -2.7717 55 3.4e-03
Maturation of SSU-rRNA -2.7465 52 3.4e-03
Mitochondrial gene expression -2.7446 54 3.4e-03
Organelle fission -2.7377 226 6.6e-03
Ribosome assembly -2.7358 77 3.4e-03
Regulation of nuclear division -2.735 79 3.4e-03
Maturation of LSU-rRNA -2.713 42 3.4e-03 
  pathwayListData.out = pathwayListData() 
 enrichmentPlot(pathwayListData.out, 25  ) 

  enrichmentNetwork(pathwayListData.out )  

  enrichmentNetworkPlotly(pathwayListData.out) 

   PGSEAplot() # pathway analysis using PGSEA 
## Error in findContrastSamples(input_selectContrast1, colnames(convertedData.out), : object 'c.out' not found

9. Chromosome

 input_selectContrast2 = limma.out$comparisons[1] 
 #input_selectContrast2 = limma.out$comparisons[3] # manually set
 input_limmaPvalViz <- 0.1  #FDR to filter genes
 input_limmaFCViz <- 2  #FDR to filter genes 
 genomePlotly() # shows fold-changes on the genome 
## Warning in eval(quote(list(...)), env): NAs introduced by coercion
## Warning in genomePlotly(): NAs introduced by coercion

10. Biclustering

 input_nGenesBiclust <- 1000    #Top genes for biclustering
 input_biclustMethod <- 'BCCC()'    #Method: 'BCCC', 'QUBIC', 'runibic' ... 
 biclustering.out = biclustering()  # run analysis

 input_selectBicluster <- 1 #select a cluster 
 biclustHeatmap()   # heatmap for selected cluster 

 input_selectGO4 <- 'GOBP'  #Gene set category 
 # Read pathway data again 
 GeneSets.out <-readGeneSets( geneSetFile,
    convertedData.out, input_selectGO4,input_selectOrg,
    c(input_minSetSize, input_maxSetSize)  )  
 results <- geneListBclustGO()  #Enrichment analysis for k-Means clusters   
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
adj.Pval Genes Pathways
2.3e-47 106 Cell wall organization or biogenesis
1.6e-37 85 External encapsulating structure organization
1.6e-37 83 Cell wall organization
7.4e-26 52 Plant-type cell wall organization or biogenesis
1.6e-24 75 Cell differentiation
9.3e-23 85 Cellular developmental process
1.9e-22 45 Cell wall biogenesis
2.6e-22 94 Carbohydrate metabolic process
9.2e-22 32 Phenylpropanoid metabolic process
5.0e-21 64 Root development

11. Co-expression network

 input_mySoftPower <- 5 #SoftPower to cutoff
 input_nGenesNetwork <- 1000    #Number of top genes
 input_minModuleSize <- 20  #Module size minimum 
 wgcna.out = wgcna()   # run WGCNA  
## Warning: executing %dopar% sequentially: no parallel backend registered
##    Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1      1    0.270 6.450          0.606     925       938    960
## 2      2    0.224 2.470          0.415     863       885    926
## 3      3    0.475 1.110          0.908     809       837    894
## 4      4    0.603 1.860          0.971     762       794    865
## 5      5    0.671 1.710          0.997     719       755    838
## 6      6    0.652 1.560          0.976     681       718    813
## 7      7    0.722 1.450          0.985     646       684    789
## 8      8    0.663 1.500          0.929     614       653    767
## 9      9    0.705 1.260          0.962     585       624    746
## 10    10    0.733 1.200          0.950     558       597    726
## 11    12    0.774 0.983          0.953     511       548    689
## 12    14    0.819 0.984          0.939     470       505    656
## 13    16    0.834 0.887          0.903     434       467    625
## 14    18    0.919 0.760          0.977     403       433    598
## 15    20    0.913 0.725          0.948     376       402    572
## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
 softPower()  # soft power curve 

  modulePlot()  # plot modules  

  listWGCNA.Modules.out = listWGCNA.Modules() #modules
 input_selectGO5 <- 'GOBP'  #Gene set category 
 # Read pathway data again 
 GeneSets.out <-readGeneSets( geneSetFile,
    convertedData.out, input_selectGO5,input_selectOrg,
    c(input_minSetSize, input_maxSetSize)  ) 
 input_selectWGCNA.Module <- 'Entire network'   #Select a module
 input_topGenesNetwork <- 10    #SoftPower to cutoff
 input_edgeThreshold <- 0.4 #Number of top genes 
 moduleNetwork()    # show network of top genes in selected module
##  softConnectivity: FYI: connecitivty of genes with less than 4 valid samples will be returned as NA.
##  ..calculating connectivities..

 input_removeRedudantSets <- TRUE   #Remove redundant gene sets 
 results <- networkModuleGO()  #Enrichment analysis of selected module
 results$adj.Pval <- format( results$adj.Pval,digits=3 )
 kable( results, row.names=FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
adj.Pval Genes Pathways
2.3e-47 106 Cell wall organization or biogenesis
1.6e-37 85 External encapsulating structure organization
1.6e-37 83 Cell wall organization
7.4e-26 52 Plant-type cell wall organization or biogenesis
1.6e-24 75 Cell differentiation
9.3e-23 85 Cellular developmental process
1.9e-22 45 Cell wall biogenesis
2.6e-22 94 Carbohydrate metabolic process
9.2e-22 32 Phenylpropanoid metabolic process
5.0e-21 64 Root development